global r4data "/.../data"
global output "/.../output"
set more off
set linesize 255

*----------------------------------------------------------------------------------------------------
* Author: Dhiren Patki
*
* This program 
* a) estimates propensity scores 
* b) estimates regression models 
* c) prepares summary/balance tables and supplemental tables
*----------------------------------------------------------------------------------------------------

*Read in stacked cohort cells 
use "$r4data/wkerpanel_m", clear

*Cell earnings are missing if all workers retire/NILF; replace with 0
replace total_earn = 0 if total_earn==.

*Age splits
local a1=25
local a2=39
local a3=40
local a4=49
local a5=50
local a6=55
local a7=56
local a8=64
local a9=65
local a10=70

*----------------------------------------------------------------------------------------------------
*		Estimate propensity scores/define common support comparison groups
*----------------------------------------------------------------------------------------------------

*No restriction on earnings (R, LFP outcome samples)

forvalues a = 1/7{

	if `a' == 1 { //25-39
	    local amin = `a1'
	    local amax = `a2'
	}
	else if `a' == 2 { //40-49
	    local amin = `a3'
	    local amax = `a4'
	}
	else if `a' == 3 { //50-55
	    local amin = `a5'
	    local amax = `a6'	
	}
	else if `a' == 4 { //50-54
	    local amin = `a5'
	    local amax = `a6'-1	
	}
	else if `a' == 5 { //56-64
	    local amin = `a7'
	    local amax = `a8'	
	}
	else if `a' == 6 { //65-70
	    local amin = `a9'
	    local amax = `a10'	
	}
	else if `a' == 7 { //25-70
	    local amin = `a1'
	    local amax = `a10'	
	}

	**Start pscore estimation step
	preserve
		local x1 "group source_at_cm5 sex_enc age_at_c no_hsd hsd some_coll coll_plus white black asian otrace hispanic av_tearn_pre log_earnpc5 tenure_cm5 firmage_cm5 av_erngwth_pre"
		local x2 "total_earn logearn R lfp ee_sep liab_active_partcp incr_liab_partcp log_liab_active_p log_incr_liab_active_p"
		local x3 "log_size lbd_emp frac_lte45 frac_46_50 frac_51_55 frac_56_60 frac_61_65 frac_66_70 frac_gte71"
		local x "`x1' `x2' `x3'"

		keep idvar `x' numpik plan_r_age treat_flag enc_first_firmid k 
		keep if inrange(age_at_c,`amin',`amax')
		keep if k<0
		keep if tenure_cm5>=2
		gen otrace2 = asian + otrace
		drop asian otrace
		rename otrace2 otrace
		
		*Coarsen plan_r_age
		bysort idvar: egen plan_r_age_av = mean(plan_r_age)
		drop plan_r_age
		replace k = k*-1

		*Redacted: drop cases with outlying plan retirement ages 

		reshape wide `x2' `x3', i(idvar) j(k)
		rename source_at_cm5 state
		rename tenure_cm5 tenure

		local x1 "i.group i.state i.sex_enc##i.age_at_c i.tenure"
		local x2 "hsd some_coll coll_plus white black otrace hispanic"
		local x3 "plan_r_age_av log_earnpc5 av_erngwth_pre total_earn* R* lfp* ee_sep* log_liab_active_p* log_incr_liab_active_p*"
		local x4 "log_size* frac_lte45* frac_46_50* frac_51_55* frac_56_60* frac_61_65* frac_66_70*"
		local x "`x1' `x2' `x3' `x4'"
		set matsize 1000

		*Estimate pscore and generate match sample indicators
		logit treat_flag `x' 
		predict pscore, pr
		psmatch2 treat_flag, pscore(pscore) common 

		*Estimate weights
		gen temp = pscore/(1-pscore) if _support==1 //p(x)/1-p(x)
		sum temp if treat_flag==0 //preferred adjustment for ATT
		gen rho_hat = r(mean)
		gen invpwt_c = temp/rho_hat if _support==1 //ATT reweighting variable
		gen invpwt_c_size = .
		replace invpwt_c_size = numpik if treat_flag==1 & _support==1 //adjust for cell size
		replace invpwt_c_size = numpik*invpwt_c if treat_flag==0 & _support==1 //adjust for cell size

		*Check overlap condition (pscore<1 in comparison group)
		tabstat pscore if _support==1 & treat_flag==0, by(treat_flag) stats(N min max mean p50)
		tabstat pscore if _support==1 & treat_flag==1, by(treat_flag) stats(N min max mean p50)

		*Re-bin time varying variables into [c-5,c-1] averages for cells within common support
		local x1 "total_earn R lfp ee_sep liab_active_partcp incr_liab_partcp lbd_emp"
		local x2 "frac_lte45 frac_46_50 frac_51_55 frac_56_60 frac_61_65 frac_66_70 frac_gte71"
		local vars "`x1' `x2'"
		foreach v of local vars {
		      egen `v'_av = rowmean(`v'1 `v'2 `v'3 `v'4 `v'5) if _support==/*redacted*/
		}
		
		forvalues v = 1/5{
		  gen log_total_earn`v' = log(total_earn`v')
		}

		local x3 "log_liab_active_p log_incr_liab_active_p log_size log_total_earn"
		foreach v of local x3{
		      egen `v'_av = rowmean(`v'1 `v'2 `v'3 `v'4 `v'5) if _support==/*redacted*/
		}
		
		*Test for balance 
		gen male = (sex_enc==2)
		local x1 "age_at_c male no_hsd hsd some_coll coll_plus white black otrace hispanic" 
		local x2 "av_tearn_pre av_erngwth_pre total_earn_av log_total_earn_av tenure R_av lfp_av ee_sep_av"
		local x3 "log_liab_active_p_av log_incr_liab_active_p_av log_size_av plan_r_age_av"
		local x4 "frac_lte45_av frac_46_50_av frac_51_55_av frac_56_60_av frac_61_65_av frac_66_70_av frac_gte71_av"
		local persv "`x1' `x2'"
		local persv_count: word count `persv'
		local firmv "`x3' `x4'"
		local firmv_count: word count `firmv'
		local totv_count = `persv_count'+`firmv_count'+1
		matrix diff_u`a' = J(`totv_count',4,.)
		matrix diff_w`a' = J(`totv_count',4,.)

		local j = 1
		foreach v of local persv {
		    reg `v' treat_flag if _support==1 [aw=numpik], vce(cluster enc_first_firmid) //unweighted difference within common support sample
		    matrix b = e(b)
		    matrix V = e(V)
		    matrix diff_u`a'[`j',1]=b[1,2] //comparison mean
		    matrix diff_u`a'[`j',2]=b[1,1] //difference in treated gp
		    matrix diff_u`a'[`j',3]=sqrt(V[1,1]) //se
		    matrix diff_u`a'[`j',4]=2*(1-normal(abs(b[1,1]/sqrt(V[1,1])))) //pval
		    reg `v' treat_flag if _support==1 [aw=invpwt_c_size], vce(cluster enc_first_firmid) //weighted difference within common support sample
		    matrix b = e(b)
		    matrix V = e(V)
		    matrix diff_w`a'[`j',1]=b[1,2] //comparison mean
		    matrix diff_w`a'[`j',2]=b[1,1] //difference in treated gp
		    matrix diff_w`a'[`j',3]=sqrt(V[1,1]) //se
		    matrix diff_w`a'[`j',4]=2*(1-normal(abs(b[1,1]/sqrt(V[1,1])))) //pval
		    local j = `j'+1
		}
		
		local j = `persv_count' + 1
		foreach v of local firmv {
		    reg `v' treat_flag if _support==1 [aw=numpik], vce(cluster enc_first_firmid) //unweighted difference within common support sample
		    matrix b = e(b)
		    matrix V = e(V)
		    matrix diff_u`a'[`j',1]=b[1,2] //comparison mean
		    matrix diff_u`a'[`j',2]=b[1,1] //difference in treated gp
		    matrix diff_u`a'[`j',3]=sqrt(V[1,1]) //se
		    matrix diff_u`a'[`j',4]=2*(1-normal(abs(b[1,1]/sqrt(V[1,1])))) //pval
		    reg `v' treat_flag if _support==1 [aw=invpwt_c_size], vce(cluster enc_first_firmid) //weighted difference within common support sample
		    matrix b = e(b)
		    matrix V = e(V)
		    matrix diff_w`a'[`j',1]=b[1,2] //comparison mean
		    matrix diff_w`a'[`j',2]=b[1,1] //difference in treated gp
		    matrix diff_w`a'[`j',3]=sqrt(V[1,1]) //se
		    matrix diff_w`a'[`j',4]=2*(1-normal(abs(b[1,1]/sqrt(V[1,1])))) //pval
		    local j = `j'+1
		}
		*Sample sizes (cell-cohorts pooled over [c-5,c-1] PIK-YEARS)
		local totv_count_m1 = `totv_count'-1
		sum treat_flag if _support==1 & treat_flag==0 [fw=numpik]
		 matrix diff_u`a'[`totv_count_m1',1]=r(N)
		 matrix diff_w`a'[`totv_count_m1',1]=r(N)
		sum treat_flag if _support==1 & treat_flag==1 [fw=numpik]
		 matrix diff_u`a'[`totv_count',1]=r(N)		 
		 matrix diff_w`a'[`totv_count',1]=r(N)

		rename invpwt_c_size invpwt_c_size`a'
		keep idvar invpwt_c_size`a'
		save $r4data/pscore_samp`a', replace
	restore
	**end pscore estimation step
}

*Positive earnings sample (logearn, total_earn, ee_sep outcome samples)
forvalues a = 1/7{

	if `a' == 1 { //25-39
	    local amin = `a1'
	    local amax = `a2'
	}
	else if `a' == 2 { //40-49
	    local amin = `a3'
	    local amax = `a4'
	}
	else if `a' == 3 { //50-55
	    local amin = `a5'
	    local amax = `a6'	
	}
	else if `a' == 4 { //50-54
	    local amin = `a5'
	    local amax = `a6'-1	
	}
	else if `a' == 5 { //56-64
	    local amin = `a7'
	    local amax = `a8'	
	}
	else if `a' == 6 { //65-70
	    local amin = `a9'
	    local amax = `a10'	
	}
	else if `a' == 7 { //25-70
	    local amin = `a1'
	    local amax = `a10'	
	}

	**Start pscore estimation step
	preserve
		local x1 "group source_at_cm5 sex_enc age_at_c no_hsd hsd some_coll coll_plus white black asian otrace hispanic av_tearn_pre log_earnpc5 tenure_cm5 firmage_cm5 av_erngwth_pre"
		local x2 "total_earn logearn R lfp ee_sep liab_active_partcp incr_liab_partcp log_liab_active_p log_incr_liab_active_p"
		local x3 "log_size lbd_emp frac_lte45 frac_46_50 frac_51_55 frac_56_60 frac_61_65 frac_66_70 frac_gte71"
		local x "`x1' `x2' `x3'"

		keep idvar `x' numpik plan_r_age treat_flag enc_first_firmid k 
		keep if inrange(age_at_c,`amin',`amax')
		keep if k<0
		keep if tenure_cm5>=2
		keep if total_earn>0 //positive earnings restriction is imposed here!
		gen otrace2 = asian + otrace
		drop asian otrace
		rename otrace2 otrace
		
		*Coarsen plan_r_age
		bysort idvar: egen plan_r_age_av = mean(plan_r_age)
		drop plan_r_age

		replace k = k*-1

		*Redacted: drop cases with outlying plan retirement ages 
	
		reshape wide `x2' `x3', i(idvar) j(k)
		rename source_at_cm5 state
		rename tenure_cm5 tenure

		local x1 "i.group i.state i.sex_enc##i.age_at_c i.tenure"
		local x2 "hsd some_coll coll_plus white black otrace hispanic"
		local x3 "plan_r_age_av log_earnpc5 av_erngwth_pre logearn* R* lfp* ee_sep* log_liab_active_p* log_incr_liab_active_p*"
		local x4 "log_size* frac_lte45* frac_46_50* frac_51_55* frac_56_60* frac_61_65* frac_66_70*"
		local x "`x1' `x2' `x3' `x4'"
		set matsize 1000

		*Estimate pscore and generate match sample indicators
		logit treat_flag `x' 
		predict pscore, pr
		psmatch2 treat_flag, pscore(pscore) common 
		*Estimate weights
		gen temp = pscore/(1-pscore) if _support==1 //p(x)/1-p(x)
		sum temp if treat_flag==0 //preferred adjustment for ATT
		gen rho_hat = r(mean)
		gen invpwt_c = temp/rho_hat if _support==1 //ATT reweighting variable
		gen invpwt_c_size = .
		replace invpwt_c_size = numpik if treat_flag==1 & _support==1 //adjust for cell size
		replace invpwt_c_size = numpik*invpwt_c if treat_flag==0 & _support==1 //adjust for cell size

		*Check overlap condition (pscore<1 in comparison group)
		tabstat pscore if _support==1 & treat_flag==0, by(treat_flag) stats(N min max mean p50)
		tabstat pscore if _support==1 & treat_flag==1, by(treat_flag) stats(N min max mean p50)

		*Re-bin time varying variables into [c-5,c-1] averages for cells within common support
		local x1 "total_earn R lfp ee_sep liab_active_partcp incr_liab_partcp lbd_emp"
		local x2 "frac_lte45 frac_46_50 frac_51_55 frac_56_60 frac_61_65 frac_66_70 frac_gte71"
		local vars "`x1' `x2'"
		foreach v of local vars {
		      egen `v'_av = rowmean(`v'1 `v'2 `v'3 `v'4 `v'5) if _support==/*redacted*/
		}
		
		forvalues v = 1/5{
		  gen log_total_earn`v' = log(total_earn`v')
		}

		local x3 "log_liab_active_p log_incr_liab_active_p log_size log_total_earn"
		foreach v of local x3{
		      egen `v'_av = rowmean(`v'1 `v'2 `v'3 `v'4 `v'5) if _support==/*redacted*/
		}

		rename invpwt_c_size invpwt_erc_size`a'
		keep idvar invpwt_erc_size`a'
		save $r4data/pscore_ersamp`a', replace
	restore
	**end pscore estimation step
}


*Merge pscores onto main file
forvalues a = 1/7{
  merge m:1 idvar using $r4data/pscore_samp`a'
  rename _merge merge_age`a'
}

forvalues a = 1/7{
  merge m:1 idvar using $r4data/pscore_ersamp`a'
  rename _merge merge_age_er`a'
}

*----------------------------------------------------------------------------------------------------
*			Retirement, LFP, earnings, and switching outcomes
*----------------------------------------------------------------------------------------------------


*Covariates
gen age2 = age^2
gen D_m5m3 = (D_m5==1 | D_m4==1 | D_m3==1)
gen D_m2m1 = (D_m2==1 | D_m1==1)
gen D_p1p2 = (D_p1==1 | D_p2==1)
gen D_p3p4 = (D_p3==1 | D_p4==1)
gen D_p5p6 = (D_p5==1 | D_p6==1)
gen D_p7p8 = (D_p7==1 | D_p8==1)
gen D_p9p10 = (D_p9==1 | D_p10==1)
gen D_p11p12 = (D_p11==1 | D_p12==1)
local x1 "D_m5 D_m4 D_m3 D_m2 D_m1 D_0 D_p1 D_p2 D_p3 D_p4 D_p5 D_p6 D_p7 D_p8 D_p9 D_p10 D_p11 D_p12 D_p13"
local x2 "age age2 i.sex_enc i.sex_enc#c.age i.sex_enc#c.age2 hsd some_coll coll_plus white black asian hispanic log_earnpc5 av_erngwth_pre"
local x "`x1' `x2'"
local x3 "D_m5m3 D_m2m1 D_0 D_p1p2 D_p3p4 D_p5p6 D_p7p8 D_p9p10 D_p11p12 D_p13"
local xcons "`x3' `x2'"
local fe "group#k source_at_cm5 tenure_cm5 enc_first_firmid"

*Depvars
local depvars1 "lfp R"
local depvars2 "logearn ee_sep" //these variables require >0 earnings


foreach y of local depvars1{

	
	**********
	* 20-39
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a1',`a2') [aw=invpwt_c_size1], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-39 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age1 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age1 [aw=invpwt_c_size1], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-39 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age1 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-39 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age1 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-39 (cell weight))
	
	**********
	* 40-49
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a3',`a4') [aw=invpwt_c_size2], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(40-49 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age2 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age2 [aw=invpwt_c_size2], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age2 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(40-49 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age2 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (cell weight))		

	**********
	* 50-55
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a5',`a6') [aw=invpwt_c_size3], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-55 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age3 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age3 [aw=invpwt_c_size3], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-55 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age3 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-55 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age3 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-55 (cell weight))
	
	**********
	* 50-54
	**********
	local b = `a6'-1
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a5',`b') [aw=invpwt_c_size4], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-54 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age4 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age4 [aw=invpwt_c_size4], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-54 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age4 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-54 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age4 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-54 (cell weight))
	
	**********
	* 56-64
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a7',`a8') [aw=invpwt_c_size5], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(56-64 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age5 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age5 [aw=invpwt_c_size5], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(56-64 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age5 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(56-64 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age5 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(56-64 (cell weight))

	**********
	* 65-70
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a9',`a10') [aw=invpwt_c_size6], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(65-70 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age6 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age6 [aw=invpwt_c_size6], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(65-70 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age6 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(65-70 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age6 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(65-70 (cell weight))

	**********
	* 25-70
	**********	
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a1',`a10') [aw=invpwt_c_size7], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-70 (IPW))
	if "`y'"== "lfp"{
	    gen pw_samp_age7 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age7 [aw=invpwt_c_size7], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-70 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_samp_age7 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-70 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_samp_age7 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-70 (cell weight))
	
}


*Depvars (log earnings and employer switching)

foreach y of local depvars2{

	if "`y'" == "logearn"{ 
	
	**********
	* 25-39
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a1',`a2') [aw=invpwt_erc_size1], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-39 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age1 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age1 [aw=invpwt_erc_size1], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-39 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age1 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-39 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age1 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-39 (cell weight))
	
	**********
	* 40-49
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a3',`a4') [aw=invpwt_erc_size2], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(40-49 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age2 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age2 [aw=invpwt_erc_size2], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age2 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(40-49 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age2 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (cell weight))		

	**********
	* 50-55
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a5',`a6') [aw=invpwt_erc_size3], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-55 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age3 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age3 [aw=invpwt_erc_size3], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-55 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age3 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-55 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age3 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (cell weight))
	
	**********
	* 50-54
	**********
	local b = `a6'-1
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a5',`b') [aw=invpwt_erc_size4], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-54 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age4 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age4 [aw=invpwt_erc_size4], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-54 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age4 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-54 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age4 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-54 (cell weight))
	
	**********
	* 56-64
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a7',`a8') [aw=invpwt_erc_size5], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(56-64 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age5 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age5 [aw=invpwt_erc_size5], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(56-64 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age5 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(56-64 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age5 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(56-64 (cell weight))

	**********
	* 65-70
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a9',`a10') [aw=invpwt_erc_size6], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(65-70 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age6 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age6 [aw=invpwt_erc_size6], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(65-70 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age6 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(65-70 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age6 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(65-70 (cell weight))

	**********
	* 25-70
	**********	
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a1',`a10') [aw=invpwt_erc_size7], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-70 (IPW))
	if "`y'"== "logearn"{
	    gen pw_ersamp_age7 = e(sample)
	}

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age7 [aw=invpwt_erc_size7], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-70 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age7 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-70 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age7 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-70 (cell weight))

	}

	else if "`y'"=="ee_sep"{

	**********
	* 25-39
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a1',`a2') & pw_ersamp_age1 [aw=invpwt_erc_size1], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-39 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age1 [aw=invpwt_erc_size1], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-39 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age1 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-39 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age1 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-39 (cell weight))
	
	**********
	* 40-49
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a3',`a4') & pw_ersamp_age2 [aw=invpwt_erc_size2], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(40-49 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age2 [aw=invpwt_erc_size2], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age2 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(40-49 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age2 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (cell weight))		

	**********
	* 50-55
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a5',`a6') & pw_ersamp_age3 [aw=invpwt_erc_size3], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-55 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age3 [aw=invpwt_erc_size3], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-55 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age3 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-55 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age3 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(40-49 (cell weight))
	
	**********
	* 50-54
	**********
	local b = `a6'-1
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a5',`b') & pw_ersamp_age4 [aw=invpwt_erc_size4], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-54 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age4 [aw=invpwt_erc_size4], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-54 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age4 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(50-54 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age4 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(50-54 (cell weight))
	
	**********
	* 56-64
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a7',`a8') & pw_ersamp_age5 [aw=invpwt_erc_size5], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(56-64 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age5 [aw=invpwt_erc_size5], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(56-64 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age5 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(56-64 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age5 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(56-64 (cell weight))

	**********
	* 65-70
	**********
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a9',`a10') & pw_ersamp_age6 [aw=invpwt_erc_size6], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(65-70 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age6 [aw=invpwt_erc_size6], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(65-70 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age6 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(65-70 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age6 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(65-70 (cell weight))

	**********
	* 25-70
	**********	
	*p-score weighted
	reghdfe `y' `x' if inrange(age_at_c,`a1',`a10') & pw_ersamp_age7 [aw=invpwt_erc_size7], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-70 (IPW))

	*p-score weighted, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age7 [aw=invpwt_erc_size7], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-70 (IPW))

	*no re-weighting
	reghdfe `y' `x' if pw_ersamp_age7 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(25-70 (cell weight))

	*no re-weighting, consolidating coefficients
	reghdfe `y' `xcons' if pw_ersamp_age7 [aw=numpik], a(`fe') vce(cluster enc_first_firmid)
	outreg2 `x3' using $output/`y'_effects_consol, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust)) ctitle(25-70 (cell weight))

	}
		
}

*----------------------------------------------------------------------------------------------------
*Save file for disclosure statistics
*----------------------------------------------------------------------------------------------------

preserve

	keep group firmid year sex_enc age_at_c source_at_cm5 tenure_cm5 idvar enc_first_firmid D* treat_flag ///
	pw_samp_age* pw_ersamp_age* lbd_emp lbd_pay numpik invpwt_c_size* invpwt_erc_size*
	save $r4data/analysis_sampflag, replace

restore

*----------------------------------------------------------------------------------------------------
*			Heterogeneity for the 56-64 age group
*----------------------------------------------------------------------------------------------------


*****************************
* Gender
*****************************

*Covariates
local x1 "D_m5 D_m4 D_m3 D_m2 D_m1 D_0 D_p1 D_p2 D_p3 D_p4 D_p5 D_p6 D_p7 D_p8 D_p9 D_p10 D_p11 D_p12 D_p13"
local x2 "age age2 hsd some_coll coll_plus white black asian hispanic log_earnpc5 av_erngwth_pre"
local x "`x1' `x2'"
local fe "group#k source_at_cm5 tenure_cm5 enc_first_firmid"

local depvars3 "lfp"
foreach y of local depvars3{

	//Men
	*Pscore weighted
	reghdfe `y' `x' if pw_samp_age5==1 & sex_enc==2 [aw=invpwt_c_size5], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects_g, replace auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addnote(Standard errors clustered at the firm level) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(Men 56-64 (IPW))

	//Women
	*Pscore weighted
	reghdfe `y' `x' if pw_samp_age5==1 & sex_enc==1 [aw=invpwt_c_size5], a(`fe') vce(cluster enc_first_firmid)
	test D_m5=D_m4=D_m3=D_m2=D_m1=0
	local F = r(F)
	local p = r(p)
	outreg2 `x1' using $output/`y'_effects_g, auto(4) ///
	addtext(Firm FE, Yes, Group-period FE, Yes, State FE, Yes, c-5 Tenure FE, Yes) ///
	addstat(Number of firms, e(N_clust), F test, `F', Prob>F, `p') ctitle(Women 56-64 (IPW))
	
}

*-------------------------------------------------------------------------------
* Summary statistics table 
*-------------------------------------------------------------------------------

local wts "w u" //pscore weighted or cell weighted

forvalues a = 1/7{ //loop over each age bin

	foreach w of local wts{ //loop over type of weight

		clear
		svmat diff_`w'`a'
		rename diff_`w'`a'1 ctrl_mean
		rename diff_`w'`a'2 diff
		rename diff_`w'`a'3 se
		rename diff_`w'`a'4 pval

		gen v = ""
		replace v = "Age at c" if _n==1
		replace v = "Male" if _n==2
		replace v = "Less than high school" if _n==3
		replace v = "High school" if _n==4
		replace v = "Some college" if _n==5
		replace v = "College or more" if _n==6
		replace v = "White" if _n==7
		replace v = "Black" if _n==8
		replace v = "Other race" if _n==9
		replace v = "Hispanic ethnicity" if _n==10
		replace v = "Earnings pre-c-5" if _n==11
		replace v = "Earnings growth pre-c-5" if _n==12
		replace v = "Earnings" if _n==13
		replace v = "Log earnings" if _n==14
		replace v = "Tenure at c-5" if _n==15
		replace v = "Retired" if _n==16
		replace v = "In labor force" if _n==17
		replace v = "Switched c-5 employer" if _n==18
		replace v = "Log pension liability (active)/active ptcp" if _n==19
		replace v = "Log pension accrual/active ptcp" if _n==20
		replace v = "Log firm size" if _n==21
		replace v = "Avg benefit claim age" if _n==22
		replace v = "Fraction workforce <=45" if _n==23
		replace v = "Fraction workforce [46,50]" if _n==24
		replace v = "Fraction workforce [51,55]" if _n==25
		replace v = "Fraction workforce [56,60]" if _n==26
		replace v = "Fraction workforce [61,65]" if _n==27
		replace v = "Fraction workforce [66,70]" if _n==28
		replace v = "Fraction workforce >=71" if _n==29
		replace v = "Control group PIK-years" if _n==30
		replace v = "Treated group PIK-years" if _n==31
		rename v variable
		order variable ctrl_mean diff se pval
		export delimited $output/wbal`w'`a'.csv, replace

	}
}


*----------------------------------------------------------------------------------------------------
*Retirement rates computed as analogs to HRS 2004 
*----------------------------------------------------------------------------------------------------

use $r4data/disc_piklevel, clear
merge m:1 group firmid year sex_enc age_at_c source_at_cm5 tenure_cm5 using $r4data/analysis_sampflag

gen age_gp=.
replace age_gp=/*redacted*/ if pw_samp_age3==/*redacted*/
replace age_gp=/*redacted*/ if pw_samp_age5==/*redacted*/
replace age_gp=/*redacted*/ if pw_samp_age6==/*redacted*/
drop if age_gp==.

keep if group==2004 & treat_flag==/*redacted*/ //2004 control group
*working at DB sponsoring firm in 2004
gen working04= (lfp==/*redacted*/ & ee_sep==/*redacted*/ & year==2004)
bysort pik: egen temp = max(working04)
replace working04 = temp if working04==/*redacted*/
keep if working04==/*redacted*/
keep if year>2004
collapse (mean) avRrate=R (semean) stderr=R, by(year age_gp)

*Reshape 
reshape wide avRrate stderr, i(year) j(age_gp)
rename avRrate1 Rav_50_55
rename stderr1 se_50_55
rename avRrate2 Rav_56_64
rename stderr2 se_56_64
rename avRrate3 Rav_65_70
rename stderr3 se_65_70
export delimited $output/lehd_R_rates04.csv, replace



